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Abstract. The purpose of this study is to explore three numerical approaches to the 
elastic homogenization of disordered masonry structures with moderate meso/macro- 
lengthscale ratio. The methods investigated include a representative of perturbation 
methods, the Karhunen-Loeve expansion technique coupled with Monte-Carlo simula- 
tions and a solver based on the Hashin-Shtrikman variational principles. In all cases, 
parameters of the underlying random field of material properties are directly derived 
from image analysis of a real- world structure. Added value as well as limitations of 
individual schemes are illustrated by a case study of an irregular masonry panel. 
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1. Introduction 

The last decade has witnessed rapid advances in modelling and simulation of masonry 
structures, mainly in connection with reconstruction and rehabilitation of historical mon- 
uments. The major impetus for these developments came from a seminal contribution 
by Anthoine [1] , which clearly demonstrated that the complex overall behavior of elastic 
regular masonry can be systematically addressed in the framework of homogenization 
theory for periodic heterogeneous media. These results were subsequently extended to- 
wards tools capable of predicting non-linear response of masonry structures with deter- 
ministic geometry and sufficiently small ratio between the characteristic size of masonry 
bond (mesoscale) and the structural level (macroscale), see [25, 28] for an extensive 
overview of the field. Both these assumptions, however, may show to be inadequate for 
historical structures, where the spatial distribution of individual constituents is random 
rather then deterministic and the typical size of a block may well become comparable to 
the macroscopic lengthscale. 

When adopting certain simplifying assumptions, a vast body of approaches is currently 
available for the treatment of irregular masonry as a random heterogeneous material. 
Under the hypothesis of widely separated lengthscales, the well-established tools of sto- 
chastic continuum micromechanics can be adopted, in which the inhomogeneous body 
in question is replaced with a homogeneous equivalent with properties determined from 
analysis of a representative volume element, locally replacing the mesoscale level, see 
monographs [5, 41] for up-to-date reviews. To the authors' best knowledge, the only 
masonry-related study available in this field was presented by Sejnoha et al. [43] in the 
framework of stochastic re-formulation of Hashin-Shtrikman variational principles due to 
Willis [44]. Additionally, the notion of a stochastic representative element can be adopted, 
based either on matching spatial statistics as originally proposed by Povirk in [32] and 



2 



M. LOMBARDO, J. ZEMAN, M. SEJNOHA, AND G. FALSONE 



subsequently applied to masonry structures in [43, 48], or deduced from the convergence 
of apparent macroscopic properties. The latter concept was proposed by Huet [17] in 
the deterministic setting, extended by Sab [34] to random media and implemented for 
historical masonry structures by Cluni and Gusella [7, 14]. 

Complementary, systems with random material properties quantified by random vari- 
ables and deterministic (or slightly perturbed) geometry of the representative volume 
can be treated employing numerical techniques of stochastic mechanics such as Monte- 
Carlo simulations [20], perturbation-based methods [22, 35] or approaches based on an 
empirical probability distribution function [36]; see also [21] for a systematic overview. 
Most generally, uncertainties in spatial distribution and material properties of individual 
phases can be jointly characterized when resorting to random field description [42]. Under 
suitable assumptions on the underlying random field, a rigorous homogenization theory 
available in [3, 19] was used to construct efficient stochastic homogenization solvers, based 
on spectral collocation methods [18] or Fourier-Galerkin approaches [46, 47]. Recently, 
these methods were extended by Xu [45] to treat heterogeneous media with small but 
finite lengthscale contrast. 

Finally, as the macro/meso scale ratio further increases, the rapidly developing tools 
of Stochastic Finite Element (SFE) methods [2, 13, 30, 40] become applicable for the 
assessment of overall response. It should be emphasized that even though the random 
field description definitely offers a more general framework than the alternative schemes, 
its major weakness is that the random field is often introduced without a clear link with 
the heterogeneous mesostructure, see [6] for a lucid discussion on the topic. Moreover, 
the application of the random field/variable paradigms to the simulation of masonry 
structure seems to be currently missing; the only related works the authors are aware 
of is a recent contribution of Spence et al. [39] dealing with mesostructure generation of 
irregular masonry walls. 

In this contribution, three numerical approaches to the determination of the overall re- 
sponse of elastic masonry with comparable macro- and meso-lengthscales are presented. 
A unifying feature is the description of mechanical properties in the form of a random 
field with the second-order statistics consistently derived from image analysis of the in- 
vestigated structure. In Section 2, this procedure is briefiy summarized following the 
exposition of Falsone and Lombardo [10]. The second level of representation involves 
the determination of basic statistics related to the response of a finite size heterogeneous 
masonry structure. In particular, the improved perturbation method is introduced first 
in Section 3, followed in Section 4 by the Monte-Carlo approach with individual realiza- 
tions of the random field generated using the Karhunen-Loeve expansion. Section 5 is 
concerned with the application of the Hashin-Shtrikman variational principles, coupled 
with the Finite Element discretization to allow for the treatment of finite-size bodies. 
In Section 6, the results obtained with the selected methods are mutually compared on 
the basis of elastic analysis of an irregular masonry panel. Finally, Section 7 introduces 
possible extensions and refinements of the studied approaches. 

In the following text, the Voigt representation of symmetric tensorial quantities is 
systematically employed, see e.g. [4]. In particular, a, a and A denote a scalar value, a 
vector or a matrix representation of a second-order tensor and a matrix representation of 
a fourth-order tensor, respectively. 
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2. Probabilistic characterization of material property via 

mesostructural statistics 

Before getting to the heart of the matter, we begin by summarizing essential termi- 
nofogy related to the theory of random fields [42]. Given a complete probability space 
{9,JF, "P} with sample space 0, cx-algebra ^ on and probability measure "P on JF, a 
scalar random field H defined on an open set f2 C M"* is a mapping 

if:0xfi^M, (1) 

such that, for every a? G f2, H{x] 9) is a random variable with respect to the triple 
{Q^T^V}. The mean of a random field is then given as 

^iH{x) = E [H{x- 9)] = [ H{x; 9) dr{9), (2) 

for any x & Q, whereas the covariance of two random fields H and G is defined by 

Rhg{x, x') = E [{H{x; 9) - finix)) {G{x'; 9) - fiaix'))] , (3) 

with the symbol Rh{x,x') = Rhh{x,x') reserved for the autocovariance, reducing to a 
variance for x = x': 

alix) = E[{Hix;9) - f,Hix)f] . (4) 

A random field H{x; 9) is said to be homogeneous if all its joint probability distribution 
functions (PDFs) remain invariant under the translation of the coordinate system, leading 
to substantial simplification of the considered statistics: 

finix) = jiH = const, cr|^(a?) = cr|f = const, Rh{x, x') = Rh{x — x'). (5) 

Assuming that the autocovariance function can be well-approximated by an exponential 
function, the correlation length Xh is defined by means of inequality: 

2 

V||x - a^'ll > Xh : Rnix - x') < ^f-, (6) 

exp(l) 

hence quantifying the characteristic dimension of the spatial fluctuations. Finally, a ran- 
dom field is ergodic if all information on joint PDFs are available from a single realization 
of the field. 

With reference to the quantification of morphology of random heterogeneous materials, 
variable 9 simply denotes a realization of random mesostructure drawn from the ensemble 
space of all admissible configurations. Of particular importance is the characteristic 
function related to the spatial distribution of the i-th phase: 

^ ^ ' ^ \ otherwise, ^ ' 

where il^^\9) is the domain occupied by the i-th. phase for realization 9 and i can take 
values {s,m}, where s denotes the stone phase and m refers to the mortar phase. The 
characteristic functions of individual phases are not independent, once e.g. the "stone" 
characteristic function is provided, the complementary descriptor follows from 

X^"^\x;9) + x^^\x;9) = l. (8) 

Therefore, we concentrate on the stone phase in the sequel. 
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When assuming the statistically uniform and ergodic media, the basic spatial statistics 
is provided by 

/XxM=cW, R^is,{x - x') = Si'\x - x') - {c^'^)\ (9) 

where c^**-* is the volume fraction of the relevant phase and 82^'^ coincides with the two- 
point probability function, defined for generic phases i,j G {s,m} as [41] 

Si''\x,x') = E[x^^\x;9)x^^\x';9)], (10) 

hence quantifying the probability of two points x and x' being located in phases i and j 
(with 5*2*^ abbreviating 82^^). 

The statistical descriptors of real mesostructures can be evaluated on the basis of a 
digitized images of the sample, leading to the discretization of the characteristic function 
in terms of an A^^^ x Ny bitmap. Replacing the point coordinate {x,y) by the pixel 
located in the i-th row and the j-th column, the characteristic function is defined by the 
discrete value Xsihj)- The estimates of one-point and two-point correlation functions, 
under the periodic boundary condition, follow from relations (see, e.g., [11]) 

^ (11) 

. Af. Ny 

S^2'\m,n) ^ ^^;^W(,,^-);^W(l + (^ + rn)%iV,,l + (j+n)%iV,), (12) 

" 1=1 ]=1 

where m and n are distances between two generic points measured in pixels and a%b 
denotes a modulo b. Note that the sums (11) and (12) can be efficiently evaluated using 
the Fast Fourier transform techniques; see e.g. [41]. To automate the acquisition of these 
functions, a software working in MATLAB was implemented [10], covering all the basic 
steps of mesostructure quantification with the data provided in the form of a color image, 
see Figure 1 for an illustration of the procedure. 

Apart from providing basic spatial statistics of random mesostructure, the phase char- 
acteristic functions allow us to directly express the matrix- valued field of material prop- 
erties in the form 

C{x; 9) = x^'\x; 9)C^'^ + x^'^Kx; 9)0^""^ = C^"^) + x^'^a;; 9){C^''^ - C^"^)), (13) 

where G^*-* and G*-'"-* are the deterministic material stiffness matrices of the two con- 
stituents. The mean and the covariance functions then follow from Equations (2) and (3): 

mM) = C^ + c^-'iClf-Cif), (14) 

Rc,M^-x') = R^^,ix-x'){clf-C^){ci^-Cl-^). (15) 

A similar representation is available when phase properties become random variables, 
see [10] for additional discussion. 

3. Improved perturbation method 

Consider a mechanical system with the randomness in material properties specified in 
terms of the random field x^'^Kx;9). In the context of finite element analysis of static 
problems, the discretized form of equilibrium equations reads [4] 

K,ix^^\x;9))u,i9) = Fj,, (16) 
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Figure 1 . Example of the program used to obtain correlation function for 
a chaotic masonry panel (San Marco d'Alunzio, Italy). 



where h is a characteristic element size, the force vector is assumed to be determin- 
istic and the global stiffness matrix is stochastic due to uncertainty in the material 
properties, which makes the nodal displacement vector Uh non-deterministic as well. 

To handle the random field in Eq. (16) computationally, an appropriate discretiza- 
tion technique in the stochastic variable has to be used. In this work, the widely used 
mid-point method is employed to represent the random field consistently with the un- 
derlying finite element mesh. Therefore, two different considerations control the size of 
an element h, cf. [2]. The first one comprises the usual "deterministic" criteria, where 
the mesh size is governed by expected stress gradients and geometry. The additional 
requirement is linked with the correlation length A^(s); the distance between two adjacent 
random variables has to be short enough to capture the essential features of the random 
field. The general recommendation is to choose 2h ~ X^(s) to describe the stochastic field 
with sufficient accuracy [29]. 

In the current implementation, the value at the element center is used to characterize 
the stochastic field, thus yielding a representation in the form of a vector of random 
variables 



with xjfg(6') being the value of x^'^\'^\(^) fhe e-th element centroid and A^'e denoting 
the number of elements. The element stiffness matrix is calculated from the standard 
finite element methodology and is expressed as [4] 



Ke (x£(^)) = / S./(^)Ce (xK(^)) Sv(^) da.. 
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where B^^e is the deterministic displacement-to-strain matrix related to the e-th element 
and the material stiffness matrix Ce follows from Eq. (13). After the assembly procedure, 
the global form of equilibrium equations becomes 

K,{x'^^\e))u,{e) = F,. (19) 

Among the various perturbative SFE approaches proposed in literature, an improved 
perturbation technique proposed by Elishakoff et al. [9] is employed in this work. When 
compared to the traditional first-order expansion schemes, the added value of the adopted 
method is that the mean value of the response variables depends on the covariance infor- 
mation on uncertain input parameters, thereby optimally utilizing the available second- 
order statistics, see also [24, Chapter 5] for further discussion. Following this approach, 
the mean of the response vector Uh is given by 

ti^^ = A~,'Fh, (20) 

where 



with the stiffness matrix sensitivities provided by 



A, = K,,o -EE ^^'Wf'^'^'^ {Kkfi-') Kl^, (21) 

i=i j=i 



In addition, the autocovariance matrix of displacements follows from 

Run = K^j^ChKo^h \ (23) 

where 

^^ = Y.Y. \f\^^K'k,f^^. (^^J ^K'r (24) 
i=i j=i 

see [9] for additional details. 

4. Karhunen-Loeve expansion 

The application of Karhunen-Loeve expansion (KLE) to stochastic boundary value 
problems has been pioneered by Ghanem and his co-workers [12, 13] and provides an 
alternative way to random field generation. The KLE can be seen as a special case 
of the orthogonal series expansion where the orthogonal functions are chosen as the 
eigenfunctions of a Fredholm integral equation of the second kind with autocovariance as 
kernel [49]. 

With reference to the mesostructure-based random fields considered in the current 
work, we start from the KLE of the characteristic function x^^^K^'i ^) in the form 



X^-'\x;e) = fi^,s,ix) + J2v^iU0)hix), (25) 

i=l 

where Aj and fi{x) are the eigenvalues (decreasing in magnitude) and eigenfunctions of 
the autocovariance R^(s) {x, x'), {C,i{9)} is a set of random variables [40]. Note that the use 
of KLE is limited to the representation of input random fields as the covariance structure 
needs to be specified a priori. Since the kernel R^(s){x,x') is bounded, symmetric and 
non-negative, it has all eigenfunctions mutually orthogonal and forming a complete set 
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spanning the function space to which belongs. Therefore, the autocovariance 

function can be decomposed into 

oo 

R^,.,{x,x') = J2^iM^)M^')^ (26) 
1=1 

with eigenfunctions fi{x) and eigenvalues Aj found as the solutions of the homogeneous 
Fredholm integral equation of the second kind 

/ R^is,{x,x')Mx')dx' = XJ,{x). (27) 
Jn 

The parameter ^i{6) in Eq. (25) corresponds to an uncorrelated standardized random 
variable expressed as 




The most important aspect of the representation (25) is that the fluctuations are decom- 
posed into a set of deterministic functions in the spatial variables separately multiplying 
purely random coefficients. 

In practical implementations, the series (25) and (26) are truncated after M terms, 
yielding the approximations 

M 

X^'\x-e) ^ fi^,s,{x) + Y.'^MO)fi{^), (29) 

1=1 

M 

R^is,{x,x') ^ ^A,/,(x)/,(a;')- (30) 

i=l 

Such spatial semi-discretization is optimal in the sense that the mean square error result- 
ing from a truncation after the M-th term is minimized [13]. 

The efficiency of KLE for simulating random fields crucially hinges on accurate eigenval- 
ues and eigenfunctions of the covariance kernel. In this paper, the Galerkin method with 
orthogonal polynomial basis functions is employed to solve the Fredholm equation (27) 
for general domains and autocovariance functions, see [24, Chapter 7] for implementation 
details. In addition, a careful convergence study of truncated KLE presented in [16] has 
demonstrated, for specific classes of stochastic fields, the dependence of the optimal value 
of M on the ratio of the characteristic domain length L to the correlation parameter X^(s) . 
For weakly correlated processes {\^(s)/L <^ 1), the higher order eigenvalues cannot be 
neglected without having a serious impact on the accuracy of the simulation. 

Such behavior is illustrated by means of Figure 2, showing the decay of eigenvalues 
of KLE with the covariance kernel determined for the masonry sample presented in Sec- 
tion 2. In addition, several associated eigenfunctions are collected in Figure 3. Clearly, 
the random field under consideration is weakly correlated as \^(s) / L ^ 10/120, see Fig- 
ure 1, and a large number of terms (M 200) is needed to capture fine features of the 
covariance, cf. Figure 3(d). 

With a KLE of the spatial autocovariance function at hand, the individual realizations 
of the heterogeneous body can be efficiently generated once an appropriate model for the 
random field x*^*^(a3; 6) is adopted. In the current study, we assume that the random field 
is Gaussian, for which the coefficients ^i{9) in Eq. (28) become independent standard 
Gaussian variables of zero mean and unit variance. It should be emphasized that the 
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Figure 3. Examples of eigenfunctions /«; (a) i = 1, (b) i = 2, (c) i = 12, 
(d) i = 24. 

assumption of Gaussian form of characteristic function is somehow questionable for binary 
heterogeneous media, where the log-normal or beta probability densities appear to be 
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more appropriate to reflect the intrinsic discreteness of tlie random field. The Gaussian 
assumption is adopted here mainly due to simplicity of the resulting simulation algorithm 
based on well-established routines, see also related works [6, 45, 46] for further discussion. 

The final step of the KLE-based solver involves the determination of the response 
statistics for a structure with material stiffness determined from Eq. (13): 

C{x; 9) ^ C^-") + l^cW + J2 ^/\^mf^ix?j (C^^^ - C(™)). (31) 

The frequently adopted framework of spectral SEE [13, 30], where the response variable 
is discretized using the polynomial chaos expansion in the stochastic coordinate, is not 
applicable in the current the high number of KLE terms results in unmanageable 

number of polynomial chaos components. Therefore, similarly to e.g. [23, 37], a direct 
Monte-Carlo approach is adopted in the present study, leading to a simulation procedure 
summarized in Figure 4. Note that in the FE analysis, the characteristic element size h 
has to verify 2h < X(s) again to correctly reproduce the autocovariance function. 



Determination of the eigenvalues 
and functions of the 
covariance kernel 

A, and f^{x) (i = 1,2,...,M) 










Generation of a random 
uncorrelated vector 

= {C.}, size M 








Realization of random field 
x(=)(a.;e,) = c(^) +E"i ^/A-^.(9,)/,(a.) 








Solution of FE model taking account 
of the spatial variation of x^'^x; 6j) 








Computation of mean and variance 
for the responses of the n samples 





iterations 



Figure 4. Overview of the Monte-Carlo simulation coupled with KLE. 
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Once the sampling phase is completed, the unbiased mean and covariance of displace- 
ment vectors are provided by 

1 

= -E^^(^^)' (32) 
i=i 

1 " 

^-'^ = ;^EM^.)K(^.))^-^A^^.KJT' (33) 

where, in accord with Figure 4, n denotes the number of simulations and 9j is used to 
denote the j-th deterministic realization. 

5. Hashin-Shtrikman variational approach 

The last approach investigated here builds on the classical Hashin-Shtrikman varia- 
tional principles for the heterogeneous media [15], extended to the stochastic setting by 
Willis [44] . The basic idea of the method is the introduction of a reference homogeneous 
body with stiffness tensor Go, employed in the analysis instead of an inhomogeneous 
realization C{x;6), recall Eq. (13). The heterogeneity of the material is compensated 
using the polarization stress t{x;6), resulting from the stress equivalence condition: 

fT{x; 9) = C{x; e)e{x; 9) = CqS^x; 6) + t{x- 6), (34) 

with cr and £ denoting the configuration-dependent stress and strain fields. The additional 
unknown follows from stationarity conditions 

[u{x;6),t{x;6)) = argminstat IIhs (i^(ic), ci;(£c); ^^) (35) 

v(x) w{x) 

where "arg min F" denotes the minimizer and "arg stat F" the stationary point of a func- 
tional F, respectively, and IIhs stands for the Hashin-Shtrikman (HS) energy functional 
provided by 

Uns{v{x),uj{xy,e) = \ f e'{v{x))Coe{v{x))dx- f v'{x)f{x)dx 

Jn Jn 

— / v'^ {x)t{x) dx + / s'^ {v{x))uj{x) dx 
Jvt Jn 

- ^ [ u;'^{x)[Cix;9) -CoY^uj{x)dx. (36) 

Jn 

In Eq. (35), v and are uj denote trial values of displacement field and polarization stresses, 
while f{x) are deterministic body forces and t{x) boundary tractions acting on Ft, re- 
spectively, cf. [26]. It can be shown that when the optimization in (35) is performed 
exactly, the stationary value of functional Hhs coincides with the actual energy stored in 
the system for realization 6. Moreover, it holds [15, 33] 

Hhs {uix; e),u;{x); 9) < Hhs {nix; 6), t{x; 6); 9) (37) 

whenever (C(x; 6) — Co) is positive-definite; when Co is chosen such that the difference 
becomes negative-definite, the inequality is reversed. 

The elementary statistics of displacements and polarizations associated with probabil- 
ity density 'P{9) follow directly from a stochastic variant of Eq. (35): 



= / (arg min stat IiYis{.v{x;e),uj{x;e);e)] dV{9). 
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Following the approach of Willis [44], the previous problem is solved approximately by 
considering the following ansatz for displacements and polarizations: 

9) = u^{x) + ui{x- 9), t{x- 9) ^ x^'\x\ 9)t^'\x) + 0)t^^\x), (39) 

where Uq is the deterministic displacement of the reference body subject to distributed 
body forces and tractions, Ui stores the configuration-dependent displacement due to the 
polarization stress expressed using a non-local operator Fq, cf. Eq. (42) bellow, 

ui{x;9) = - I roix,x')T{x';9)dx', (40) 
Jn 

and T*^*-* denotes the deterministic value of the polarization stress related to the i-th. phase. 

In accord with the standard Galerkin procedure, an identical form is adopted for the 
trial values of displacements and polarizations. Upon exchanging the order of optimiza- 
tion, the optimality conditions of problem (38) reduce to identity [26, 27] 

[ Lj^'^^{x)c^'^[C^'^ -Co]~\^'\x)dx + (41) 

V / Lj^'^^ix) [ si''\x-x')To{x,x')T^^\x')dx' = [ u;^'\xyc^'^e{uo{x))dx, 
■ Jn Jn Jn 

to be satisfied for an arbitrary ljW and i,j G {s,m}. 

Two levels of approximation are generally needed to fully discretize the system (41). 
The first step involves discretizing the Fq operator together with the "reference" strain 
distribution, which in the context of the adopted Finite Element approximation be- 
come [26, 27] 

Toix, x') ^ Bl{x)K^lBf{x'), e{uo{x)) ^ Bl{x)uo,h, Uo{x) ^ NI{x)uom, (42) 

where, in analogy with Section 3, Khfl denotes the stiffness matrix of the reference 
structure, N"^ is the matrix of shape functions and Bh = dN"^ is the displacement-to- 
strain matrix and ito,h stands for nodal displacement vector determined for the reference 
problem [4]. In the second step, phase polarization fields are parameterized in the form, 
cf. [26] 

r^H^)^Nl{x)4\ (43) 

where A/"^ is the matrix of shape functions to approximate the polarization stresses. It is 
worth mentioning that a detailed one-dimensional study presented in [38] demonstrated 
that, similarly to the remaining approaches, the characteristic element size 2h ^ X^(s) is 
again necessary to achieve a sufficient accuracy of the obtained statistics. 

Employing the approximations (42) and (43), the stationary conditions (41) yield the 
system of hnear equations 

K«d« + $:K^^-^4^-) = il«, (44) 
j 
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with the individual terms provided by 

Kf = ! NI\x)c^^[C^^-C,]-'nI{x)Ax, (45) 

= / I Nl\x)S^^'\x~x')T,,u{x,x')Nl{x')dxdx\ (46) 
Jn Jn 

= [ NjJ{x)c^^^Bl{x)uo,,dx. (47) 
Jn 

Once the degrees of freedom related to the phase polarization stresses are determined 
from system (44), the mean of displacement value becomes [26] 

^^^^ix) = Nl{x) (^,^,,-K,xj^Bl\x')ti^^{x')dx^ , (48) 

with 

M.,(x') = Nl{x') (c(-)4") + cWdi^)) . (49) 

In addition to the mean response, the HS approach offers an alternative way to estab- 
lishing confidence-like bounds on the expected displacements by varying the auxiliary 
stiffness Cq. In particular, it follows from Eq. (37) that selecting the reference medium 
such that Co = minj(C*-*'') yields an upper bound of the stored energy (and therefore 
the upper "energetic" bounds of the displacements), whereas the choice Cq = maXj(C*-*^) 
results in a lower bound on displacements. Finally, selecting Cq such that the difference 
{C — Co) becomes indefinite provides general variational estimates of the basic statistics, 
cf. [8]. 

6. Numerical example 

In this Section, the essential features of the proposed numerical methods are illustrated 
by studying elastic response of an irregular masonry structure with dimensions shown in 
Figure 5 and constant thickness of 0.12 m. The plane stress assumptions were adopted 
in the analysis; the structure was subject to a uniform pressure applied at the top edge 
and to the self-weight (a deterministic specific gravity equal to 20 kNm~^ was assumed 
for simplicity). Material constants of individual constituents were considered to be deter- 
ministic, the concrete values of the Young moduh E^™) = 1, 200 MPa, E^'^ = 12, 500 MPa 
and of the Poisson ratios z/^™^ = 0.3 and z/*-*^ = 0.2 were selected following [7]. The ge- 
ometrical uncertainty due to irregular configuration of individual phases was quantified 
on the basis of image analysis data presented in Figure 1. 

The finite element model of the example problem was based on a regular discretization 
of the domain using 24 x 24 square bilinear elements with four integration points. Note 
that such a resolution corresponds to the element edge approximately equal to a half of 
the geometrical correlation length, which is fully consistent with general rules discussed 
in Sections 3-5. 

The results presented for the KLE-based solver were derived from n = 1,000 simu- 
lations. For simplicity, only the Young modulus E considered in the form of a random 
field (see Figure 6 for an illustration), whereas the Poisson ratio was set to a deterministic 
value determined from the Voigt estimate of the homogenized stiffness matrix 

^^ = c(")c('") + c(^)c(^). (50) 
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\ L=1.20m ^ 

Figure 5. Scheme of an illustrative example 

Finally, based on a systematic one-dimensional study of the HS approach presented in [38], 
the element-wise constant discretization of the phase polarization stresses with four-point 
integration scheme was adopted to ensure sufficient resolution of the spatial statistics. 
Additional implementation-related details can be found in [24]. 





100 120 



(a) (b) 

Figure 6. Realization of Gaussian random field of Young's modulus 

Before presenting the comparison of individual approaches, we concentrate first on 
the effect of reference media on the HS-based predictions. To this end, the expected 
values of nodal displacements are plotted in Figure 7 for several representative choices of 
Cq. In particular, owing to the dominant fraction of the stiffer phase (c*-'^-' = 67%) in the 
considered structure, the lower energetic bound can be expected to be substantially closer 
to the "true" statistics than the corresponding upper bound, which in the current case 
seems to be too inaccurate for practical use. Additional estimates can be generated by 
the Voigt-type choice (50) or by setting the reference medium to the arithmetic average 
of properties of individual constituents, the value commonly adopted in the polarization- 
based numerical method due to Mouhnec and Suquet [31]. As expected, the response 
corresponding to such choices is comparable to the lower bound and will be used in the 
sequel for the comparison with the candidate approaches. 
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Figure 7. Expected value of nodal displacements via HS-based FE analysis 

The basic statistics of nodal displacements, as predicted by different methods, are mu- 
tually compared in Figure 8. In addition, we present the confidence bounds in the form 
fi^^{x) ±(Tu^{x), determined on the basis of the second-order statistics for the improved 
perturbation method (24) or KLE (33). In general, it can be seen that the perturba- 
tive method leads to a substantially wider confidence interval when compared to the 
Monte-Carlo simulation approach, in spite of a moderate number of simulations used by 
KLE solver to estimate the overall statistics. For both methods, the confidence intervals 
remain bounded from above by the corresponding HS value. Moreover, for appropriate 
choices of the reference stiffness matrix, the HS method recovers the predictions provided 
by the alternative approaches. For the current setting, selecting Cq according to the rule 
of mixtures yields the displacement values almost identical to that of the KLE solver, 
whereas the response related to the arithmetic average well approximates the improved 
perturbation result. These results provide just another highlight of the importance of a 
proper choice of the reference media in the HS-based schemes, see e.g. [8, 38] for further 
discussion. 

The final comment concerns the computational complexity of individual approaches. It 
can be stated that the cost of the improved perturbation method and the HS-based solver 
is roughly comparable, whereas the KLE approach leads to an approximately three-fold 
increase in the simulation time. Higher cost of the latter method can be attributed to a 
large number of terms appearing in the expansion (29); the computational cost, however, 
is compensated by generality of the Monte- Carlo framework and can be further reduced 
by parallelization of the problem. 



In this contribution, the applicability of three distinct approaches to mesostructure- 
based random field simulation of irregular historic masonry was investigated. The nu- 
merical results obtained for a finite-size elastic panel allow us to reach the following 
conclusions: 

• The elements of quantification of random spatial statistics can be efficiently used 
to construct realistic first- and second-order statistics of stationary random fields. 



7. Conclusions and future work 
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Figure 8. Basic statistics of nodal displacements on the top of the panel 



• The improved perturbation method utilizes the second-order statistics when de- 
termining the mean response of the system, which generally leads to narrower 
estimates when compared with the basic method, e.g. [24, Chapter 5]. In the 
current case, however, the uncertainty in the obtained statistics is higher than for 
the KLE algorithm, mainly due to a relatively high contrast of phase stiffnesses. 

• The Karhunen-Loeve series representation coupled with the Monte Carlo approach 
provides an interesting alternative to the perturbation-based method, even at in- 
creased computational cost. When applied to realistic structures, however, a large 
number of terms seems to be necessary to capture the available covariance infor- 
mation. Moreover, the validity of the Gaussian assumption needs to be critically 
assessed. 

• The Hashin-Shtrikman approach takes advantage of the specific form of the ran- 
dom field and therefore optimally utilizes the available information. The overall 
response is in this case, however, highly dependent on the choice of the reference 
medium, for which there is no general rule. 

Even though the results of this pilot study have provided valuable insights into the pros 
and cons of individual methods, they do not allow for directly quantifying the accuracy 
of individual methods as the reference solution is not available. Similarly to a recent 
study [38] , such a comparison can be based on a synthetic mesostructural model such the 
one proposed by Spence at al. [39]. This particular topic enjoys our current interest and 
will be reported separately. 
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